function hfile2 = vmr_HiResRescale(hfile, res, imeth, cutoff)
% VMR::HiResRescale  - upsample ISO-Voxel VMR to higher resolution
%
% FORMAT:       hresvmr = vmr.HiResRescale([res, imeth])
%
% Input fields:
%
%       res         resolution (default: 0.5mm)
%       imeth       interpolation: {'linear'}, 'cubic', 'lanczos3'
%       cutoff      if given and false, do not cutoff all-zero planes
%
% Output fields:
%
%       hresvmr     HiRes iso-voxelated VMR
%
% Note: the framing cube will be chosen automatically:
%        512 for  0.5 ... <1.0
%       1024 for >0.0 ... <0.5
% Also: the function only supports VMR that are ISO-voxelated

% Version:  v0.7f
% Build:    9030323
% Date:     Mar-03 2009, 11:08 PM CET
% Author:   Jochen Weber, SCAN Unit, Columbia University, NYC, NY, USA
% URL/Info: http://wiki.brainvoyager.com/BVQXtools

% argument check
if nargin < 1 || ...
    numel(hfile) ~= 1 || ...
   ~isBVQXfile(hfile, 'vmr')
    error( ...
        'BVQXfile:BadArgument', ...
        'Invalid call to ''%s''.', ...
        mfilename ...
    );
end
bc = bvqxfile_getcont(hfile.L);
if ~isfield(bc, 'VoxResX') || ...
    any([bc.VoxResX, bc.VoxResY, bc.VoxResZ] ~= 1) || ...
    bc.FramingCube ~= 256
    error( ...
        'BVQXfile:InvalidObject', ...
        'Only defined for 1mm ISOvoxel VMRs in a 256 FramingCube.' ...
    );
end
if nargin < 2 || ...
   ~isa(res, 'double') || ...
    numel(res) ~= 1 || ...
    isinf(res) || ...
    isnan(res) || ...
    res < 0.1 || ...
    res >= 1
    res = 0.5;
end
if nargin < 3 || ...
   ~ischar(imeth) || ...
    isempty(imeth) || ...
   ~any(strcmpi(imeth(:)', {'linear', 'spline', 'cubic', 'lanczos3'}))
    imeth = 'linear';
else
    imeth = lower(imeth(:)');
end
if res < 0.5
    fcube = 1024;
else
    fcube = 512;
end
fch = fcube / 2;
fcm = fch - 1;
if nargin < 4 || ...
   ~islogical(cutoff) || ...
    isempty(cutoff) || ...
    cutoff(1)
    cutoff = true;
else
    cutoff = false;
end

% which data
if ~isempty(bc.VMRData16) && ...
    numel(size(bc.VMRData16)) == numel(size(bc.VMRData)) && ...
    all(size(bc.VMRData16) == size(bc.VMRData))
    sdata = bc.VMRData16(:, :, :);
else
    sdata = bc.VMRData(:, :, :);
end

% find source offsets
mxx = max(max(sdata, [], 2), [], 3);
mxy = max(max(sdata, [], 1), [], 3);
mxz = max(max(sdata, [], 1), [], 2);
sox = find(mxx > 0);
soy = find(mxy > 0);
soz = find(mxz > 0);

% copy VMR, but set some fields to better values...
ndata = sdata(1);
sdata = single(sdata);
ndata(1) = 0;
hfile2 = aft_CopyObject(hfile);
sfile2 = bvqxfile_getscont(hfile2.L);
sfile2.F = '';
sfile2.C.FramingCube = fcube;
sfile2.C.VoxResX = res;
sfile2.C.VoxResY = res;
sfile2.C.VoxResZ = res;
if isempty(sox)
    if cutoff
        sfile2.C.VMRData = ndata;
        sfile2.C.VMRData(1:4, 1:4, 1:4) = 0;
        sfile2.C.OffsetX = fch - 2;
        sfile2.C.OffsetY = fch - 2;
        sfile2.C.OffsetZ = fch - 2;
        sfile2.C.DimX = 4;
        sfile2.C.DimY = 4;
        sfile2.C.DimZ = 4;
    else
        sfile2.C.VMRData = ndata;
        sfile2.C.VMRData(1:fcube, 1:fcube, 1:fcube) = 0;
        sfile2.C.OffsetX = 0;
        sfile2.C.OffsetY = 0;
        sfile2.C.OffsetZ = 0;
        sfile2.C.DimX = fcube;
        sfile2.C.DimY = fcube;
        sfile2.C.DimZ = fcube;
    end
    bvqxfile_setscont(hfile2.L, sfile2);
    return;
end

% determine source center coordinate
ccsrc = (bc.FramingCube - 1) / 2 - [bc.OffsetX, bc.OffsetY, bc.OffsetZ];

% generate sampling grid
xsmp = res / bc.VoxResX;
ysmp = res / bc.VoxResY;
zsmp = res / bc.VoxResZ;
xs = xsmp * (-0.5:(fcube - 1)) + (ccsrc(1) - fcm * xsmp) + 1;
ys = ysmp * (-0.5:(fcube - 1)) + (ccsrc(2) - fcm * ysmp) + 1;
zs = zsmp * (-0.5:(fcube - 1)) + (ccsrc(3) - fcm * zsmp) + 1;
xg = find(xs > (sox(1) - 3) & xs < (sox(end) + 3));
yg = find(ys > (soy(1) - 3) & ys < (soy(end) + 3));
zg = find(zs > (soz(1) - 3) & zs < (soz(end) + 3));

% generate output
if cutoff
    ndata(1:numel(xg), 1:numel(yg), 1:numel(zg)) = ndata(1);
else
    ndata(1:fcube, 1:fcube, 1:fcube) = ndata(1);
end
xs = [Inf; xs(xg(1)); xs(2) - xs(1); xs(xg(end))];
ys = [Inf; ys(yg(1)); ys(2) - ys(1); ys(yg(end))];
r = [xs, ys, ys];

% determine progress bar capabilities
try
    pbar = BVQXprogress;
    BVQXprogress(pbar, 'setposition', [80, 200, 640, 36]);
    BVQXprogress(pbar, 'settitle', sprintf('%s interpolating VMR...', imeth));
    BVQXprogress(pbar, 0, 'Building interpolation grid...', 'visible', 0, 1);
catch
    pbar = [];
end

% iterate over Z grid
pbv = 0;
if ~isempty(pbar)
    BVQXprogress(pbar, 0, 'Interpolating...', 'visible', 0, size(ndata, 3));
end
for zc = zg
    r([2, 4], 3) = zs(zc);
    pbv = pbv + 1;
    ndata(:, :, pbv) = flexinterpn_method(sdata, r, imeth);
    if ~isempty(pbar)
        BVQXprogress(pbar, pbv);
    end
end
if ~isempty(pbar)
    closebar(pbar);
end

% set into new VMR
if cutoff
    sfile2.C.OffsetX = round(r(2, 1) / r(3, 1));
    sfile2.C.OffsetY = round(r(2, 2) / r(3, 2));
    sfile2.C.OffsetZ = round(zs(1) / r(3, 3));
else
    sfile2.C.OffsetX = 2 * bc.OffsetX;
    sfile2.C.OffsetY = 2 * bc.OffsetY;
    sfile2.C.OffsetZ = 2 * bc.OffsetZ;
end
sfile2.C.DimX = size(ndata, 1);
sfile2.C.DimY = size(ndata, 2);
sfile2.C.DimZ = size(ndata, 3);

% simple if only 8-bit data
if ~sfile2.C.VMR8bit || ...
    isa(ndata, 'uint8')
    sfile2.C.VMRData = ndata;
    bvqxfile_setscont(hfile2.L, sfile2);
    return;
else
    sfile2.C.VMRData16 = ndata;
end

% otherwise VMR data has to be reconstructed
nh = cumsum(hist(single(ndata(:)), single(0:max(ndata(:)))));
ht = find(nh > 0.999 * numel(ndata));
if isempty(ht)
    ht = numel(nh);
end
ndata(ndata > ht(1)) = ht(1);
sfile2.C.VMRData = uint8(225 * (single(ndata) ./ single(max(ndata(:)))));
bvqxfile_setscont(hfile2.L, sfile2);
